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A radiative transfer scheme for cosmological reionization 
based on a local Eddington tensor. 
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ABSTRACT 

A radiative transfer scheme is presented, based on a moment description of the equa- 
tion of radiative transfer and the so-called "Ml closure model" for the Eddington ten- 
sor. This model features a strictly hyperbolic transport step for radiation : it has been 
implemented using standard Godunov-like techniques in a new code called ATON. 
Coupled to simple models of ionization chemistry and photo-heating, ATON is able 
to reproduce the results of other schemes on a various set of standard tests such as 
the expansion of a HII region, the shielding of the radiation by dense clumps and 
cosmological ionization by multiple sources. Being simple yet robust, such a scheme is 
intended to be naturally and easily included in grid-based cosmological fluid solvers. 

Key words: Radiative Transfer - Cosmology - Methods: Numerical, N Body. 



1 INTRODUCTION 

During the early stages of the Universe, the process of struc- 
ture formation leads to the formation of the first stars at 
a redshift z ~ 10 — 20. These primordial, metal-poor stars 
should emit a substantial amount of radiation that would be 
able to ionize the neutral gas that fills the surrounding space 
(see e.g. Barkana & Loeb |2001 ) for an extensive review). 
Afterward, it would become transparent to radiation at high 
energies, where the transition is commonly pictured as ion- 
ized bubbles that expand and percolate around sources. The 
investigation of the the distribution of neutral gas at high 
redshift is therefore a great source of knowledge on the first 
luminous objects, on the physical conditions in which they 
appear and on the cosmological context that lead to them. 
For this purpose, several experiments were or are about to 
be set up : one can mention LOFAR or SKA which aim at 
detecting the redshifted 21 cm signal which would come up 
from the neutral hydrogen. 

From a theoretical point-of-view, the interest of pursu- 
ing so-called 'full-physics' cosmological simulations has been 
recently emphasized by different groups These numerical ex- 
periments can predict the history of star formation, the large 
scale gas distribution hosting small scale disc-like objects, 
the creation and ejection of metals, etc... While intensively 
tested at low redshift, only a few observational comparisons 
are available for 2 > 3 and basically none for z > 5. In 
this context, the comparison of simulations' prediction to 



aubert@astro.u-strasbg.fr 



the observed 21 cm emission would open a new range of 
epochs during which the understanding of the formation of 
structures could be tested. 

Common features of hydrodynamical simulations are 
gravity, hydrodynamics and star formation. They are mod- 
eled self-consistently and coupled to each other, up to a cer- 
tain extent. On the other hand, the radiation is often consid- 
ered as an external homogeneous background and its impact 
on the physics is not exactly consistent with the distribution 
of sources inside the simulated box. A common approach 
consists in taking in account a diffuse radiation background 
which evolves along time according to a model of cosmolog- 
ical radiation sources ( see e.g. |Madau et al] ([1999)). The 
geometrical distribution of the sources and the propagation 
of the emitted radiation in the neutral hydrogen gas cannot 
be eluded if one is interested in the transition processes that 
occur at the reionization. Note that still at lower redshift, al- 
though the universe is fully ionized, the radiation field is still 



highly inhomogeneous down to redshift 1 ( see e.g. Madau 
[etar]|T999l )). 

A great effort has been put into developing tools that 
simulate the emission and the propagation of radiation 
in cosmological boxes and its impact on the gas physics 
through heating/cooling processes and ionization (see e.g. 
Gnedin fc Abel] ( |2001[ ), [Mellema et"al] (|2006[)). A perfect 
illustration of this interest in given by |Iliev et al.| ( |2006[ ) 
where a large number of radiative transfer codes were gath- 
ered in order to be tested on the same set of numerical ex- 
periments. Several methods (ray shooting, grid-based codes, 
Monte-Carlo) and several types of implementation (post- 
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processing, self-consistent simulations) are compared and 
this project demonstrates that different methods can achieve 
very similar results even though they greatly differ in their 
conception. We would like to stress here that in this compar- 
ison paper, only one moment-based code was used, namely 



the OTVET code iGnedin fc Abel 20011, and its use was 



only restricted to strictly periodic problems. 

In this paper, we present a new moment-based method 
to perform calculations of radiative transfer. It relies on a 
rather standard hyperbolic grid-based solver and is simple 
enough to be quickly implemented. This method relies on 
a momentum description of the transfer equation via the 
conservation of radiative energy and flux. [Gnedin fc Abel] 
( 2001 1 relied on the same type of description of the equa- 



tion of radiative transfer, where the Eddington's tensor is 
constrained by the sources' geometry, assuming an optically 
thin regime. We use here to compute the Eddington tensor 
the so-called Ml closure relation, which provides a variable 
Eddington tensor that depends only on the local radiation 
flux and intensity. Non-equilibrium ionization is also taken 
in account by solving the coupled set of equations for out- 
of-equilibrium hydrogen chemistry. The current implemen- 
tation performs only post-processing of existing simulation, 
even though it can in principle be easily coupled to a grid- 
based hydrodynamical code. 

The paper is organized as follows: the first section 
briefly introduces the model. The numerical implementa- 
tion is then presented in more details. The third section is 



devoted to the same tests performed in Iliev et al. (20061, 



on rather academic situations but also on realistic cosmo- 
logical fields. The range of applications of this new scheme 
is discussed in the last section. 



2 RADIATIVE TRANSFER AS AN 
HYPERBOLIC SYSTEM OF 
CONSERVATION LAWS WITH SOURCE 
TERMS 

The current scheme is based on a momentum description of 
the radiative transfer equation. The hierarchy of equations 
is truncated at the second order and the closure relation is 
provided by the Ml relation. 



2.1 Moments of the transfer equation 

The radiative transfer equation is given by: 

c dt 

where /i/(x, n, t) is the radiation specific intensity, k^(x, n, t) 
the absorption coefficient and »?,y(x, n, t) the source function. 
They all depend on position, angle, frequency and time. The 
absorption coefficient, in the context of ionizing radiation, 
is computed from 



(1) 



(2) 



where ct^ is the photoionization cross section, and tiho the 
neutral hydrogen density. By taking the first two momenta 
of Eq. [l] two coupled equations can be obtained: 

FIE 

— -^-HVF, = -K.cS. +S'„ (3) 
at 



at 



(4) 



These four equations set the conservation of the radiative en- 
ergy El,, the zero-th order momentum of the intensity, and 
the conservation of the radiative flux Fj^, the flrst order mo- 
mentum. The lower dimensionality of these equations, plus 
their conservative form, make them more suited to a nu- 
merical treatment. However, an expression for the pressure 
tensor (i.e. the second order momentum of the intensity) 
must be provided in order to close the system described by 
Eqs.[3]and|4] This issue is addressed in Sec. 12.41 

From now on, Eq. (|3| and Eq. Q can be modifled in 
a form more suitable to the subsequent calculations. First, 
the energy and flux densities can be replaced by number 
densities : this is easily achieved by dividing Eq. ([sjl and 
Eq. Q by a single photon energy, hv. They become : 

dt 
9F, 



- VF, 



dt 



+ c VP, 



(5) 
(6) 



where N,, is the photon number density. For sake of sim- 
plicity, we have used the same notation for the photon flux 
(resp. photon pressure tensor) than for the energy flux (resp. 
energy pressure tensor), although they differ by the factor 
hv. The source term is further divided into 2 contributions 



Si, = K + N^" where K"" = n,n^ k„{T) 



(7) 



where the flrst term is the radiation coming from stars or 
quasars and the second term is the diffuse radiation due to 
recombination from H-|- . Both radiation sources are assumed 
to be isotropic, so that no source term appears in the flux 
equation. In the test section, we will compare our method to 
ray-tracing schemes developed in the context of cosmologi- 
cal reionization ( Iliev et al.||2006 l, for which recombination 
radiation is emitted along each ray, so that recombination 
radiation is not isotropic anymore. In order to mimic the ef- 
fect of ray-tracing, we optionally solve a modified equation 
for the flux 



at 



IS u 



(8) 



We call this approximation the "ray-tracing scheme" , al- 
though the underlying method still makes use of the 2 first 
moments of the radiative transfer equation. 



2.2 Single group radiative transfer 

In the current implementation of ATON, we have restricted 
ourselves to the ionization of a single specie, namely hydro- 
gen, and discard completely the fate of helium and other ele- 
ments. It is relatively straightforward to extend our scheme 
to a more realistic chemical composition, using for exam- 
ple the multiple-frequency approach described in [Gnedin fc| 
|Abel| ( |2001| |. This is beyond the scope of this paper. We 
further simplify the problem by considering only one pho- 
ton group, namely all photons with energy greater than the 
threshold energy for hydrogen. We use throughout this pa- 



per the notations introduced by Katz et al. ( 1996 1. The num- 
ber density of hydrogen nuclei is rin = pX/rUp (for which 
we adopt X — 0.76), while the number density for neutral 
hydrogen (resp. ionized hydrogen) is noted riHo (resp. nH_|_). 
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We define the ionizing photons number density as 



N^Au. 



(9) 



Integrating Eq. ^ and Eq. ([6| over photon frequency and 
dropping the subscript Ho since we have only one photon 
group in this paper, we get 



1h 



+ VF^ = 



rec 



dt + ^ - 



-nHoCdiviV^ + iV* +iV. 



(10) 
(11) 



where we define two frequency-averaged cross-sections 
given by 

/OO /"OO 

(JuNudu and crpF-y — I avT'^Av. (12) 

In order to simplify further the problem, we assume a simple 
reference radiation intensity, noted Joiy), for which we pre- 
compute the average cross-section as follows 



Of — UN 



a„ — r^^di/ / / — -^^di^, (13) 

"Ho "'"Ho 



where the hydrogen photoionization cross-section auf, (i^) is 



taken from Hui & Gnedin (19971. Except in section 4.1 



we will consider 10 K black-body models where aj = 
1.63 X 10~^* cm^. Likewise, the recombination radiation, 
integrated over frequency in our single energy group writes 



onH+ eH_,_ {i^, T) = nonH+ {a a — olb ) 



(14) 



where aA{T) (resp. as(r)) is the case A (resp. case B) 
recombination coefficient for H+. They are both taken from 



from Hui & Gnedin ( 19971. The set of equations we solve in 



this paper is finally 



dt 



+ VF^ = — nHoCo-^A^7 + nonH+(aA — as) 



-nHoCO"-,F^ 



(15) 
(16) 



Let us recall that this set of equation is obtained for one 
single group of photons, the ionizing ones. The same pro- 
cedure can be applied in principle to an arbitrary number 
of groups, in order to achieve a better spectral description 
of the problem. The number of systems to be solved would 
therefore scale with the number of groups considered. 



2.3 Hydrogen thermochemistry 

In order to close the last system of equation, we need to 
solve for the time evolution of the Hydrogen ionization frac- 
tion and of the gas temperature. The chemical evolution of 
neutral hydrogen is governed by a delicate balance between 
coUisional ionization, photoionization and coUisional recom- 
bination. These processes are part of the following evolution 
equation for riHo 



D 
Dt 



(«Ho) = QAnoriH^ — ISricn-aa — r^HoWHo, 



(17) 



together with charge conservation tIc — nn^ and Hydrogen 
nuclei conservation tih^ + riHo = J^h- L^Hq is the Hydrogen 
atom photoionization rate, given by (using the same nota- 
tions as for Eq. 12 1 



- 7H0 



(18) 



Radiative cooling and photoionization heating are also self- 
consistently taken into account by solving the gas thermal 
energy equation 



with e = -ritot ksT. 



(19) 



The cooling rate, £,, in erg/s/cc, is computed using standard 
coUisional cooling processes due to case A and B recom- 
bination of Hydrogen, coUisional ionization and excitation 
of Hydrogen and Bremsstrahlung. We use the cooling rates 
given by Hui & Gnedin ( 1997) , MaseUi et al.| ( |2003| ) and ref- 
erences therein. The photoionization heating rate, TL, also 
in erg/s/cc, is given by H = nugiiio where 



EHo = c 



[hi' — hi'-aQ)ui,Ni,diy' = ce^/a^/N^ 



(20) 



Following the approach used in the last section, we approx- 
imate this photoionization energy using the fiducial 10^ K 
black body radiation spectrum, so we can precompute the 
average photon energy (equals to 29.65 eV in this case). 



2.4 The Ml closure relation 

As already mentioned, we need to specify the form of the Ed- 
dington tensor in order to close the moment hierarchy and 



solve the previous set of equations. Gnedin & Abel (20011 



suggested to compute the Eddington tensor assuming an 
optically thin medium and summing up the contribution of 
all the background sources. In this way, the radiative inten- 
sity geometry was fully specified and the moment equations 
could be solved. This techniques, refereed to as the "Opti- 
cally Thin Variable Eddington Tensor" method, required to 
solve four different Poisson-like equations. Since all cosmol- 
ogy codes already have a Poisson solver to compute the dark 
matter and gas dynamics, this scheme turned out to be quite 
efficient and accurate for cosmological applications. In this 
paper, we propose to apply a very simple closure relation to 
cosmological reionization, called the "Ml approximation", 
introduced more than two decades ago to solve the radiative 
transfer equations in the optically thick limit, while retain- 



ing some accuracy in the optically thin regime (Levermore 



19841 



It relies on the assumption that the radiation angular 
distribution is axysymetric around the flux vector F, so that 
the Eddington tensor, defined as P = DA'', can be written 
in the general form: 



D 



1 



X 



I 



3x-l 



(21) 



2 2 ' 

where u is a unit vector aligned with the flux direction. We 
need also th define the reduced fiux f by: 

f-^ = /- (22) 

The Eddington factor x is a yet unknown scalar quantity 
that depends only on /, and that should satisfy 1/3 ^ x ^ 1- 
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We now need a model to specify the functional form for x(/)- 



In his review paper, Levermore ( 1984 1 discussed a great va- 



riety of closure relations. The most simple and physically 
meaningful one is called the Ml model, for which we have: 



X = 



3 + 4lfl- 



5 + 2^4-3|fl- 



(23) 



This closure relation corresponds to the angular distribution 
of a Lorentz boosted isotropic distribution ( Levermorej 1984 1 , 
as for the Cosmic Microwave Background dipole, for which 
the boost direction is ahgned with the flux vector. It has 



been shown by Dubroca & Feugeas ( 1999, ) that this closure 
relation is the only one that minimize the radiative entropy. 

This closure relation is of course a very crude approx- 
imation of the true radiation distribution. In presence, for 
example, of two distinct radiation source, the Ml model will 
replace the two sources by one "average" source in between. 
Nevertheless, this colure relation has good properties that 
we now discuss in more details. This model satisfies the phys- 
ical constraint 1/3 ^ x ^ 1- The tensor D, which describes 
the radiation's local geometry, consists in two separate con- 
tributions. The first one is an isotropic component, where 
the radiation affects all the directions in a similar way. In 
Eq. |21[ the second tensor component of D exhibits prin- 
cipal directions that are aligned with the local fiux, consis- 
tently with a free-streaming radiation. On the one hand, the 
isotropic component disappears in a pure transport regime 
with / = 1, i.e. X = 1- On the other hand, the diffusion 
regime implies / = and x ~ 1/3- Using this value in Eq. 
[21] shows that only the isotropic component remains, as ex- 
pected. In a general fashion, all the intermediate regimes 
represent local geometries where both an isotropic radiation 
and a free-streaming radiation contribute. These two limit- 
ing cases are exactly described by the Ml model, while the 
general regimes are approximated by a linear combination 
of these two limiting cases. 

The other interesting property is that this model is 
purely local, so that no expensive Poisson solvers are nec- 
essary. Even more interestingly, it can be shown that the 
left-hand side of Equations ( |15[ ) and ( |16| l defines an hyper- 
bolic system of conservation laws, with real eigenvalues cor- 
responding to waves traveling at (or close to) the speed of 
light. ( [Dubroca fc Fcugoas 1999 , Gonzalez et al.||2007[ ). We 
can therefore use standard numerical techniques designed in 



the general framework of hyperbolic conservation laws ( Toro 
|1999[ ) and apply them to cosmological reionization. 



3 NUMERICAL IMPLEMENTATION 

We now describe in details the numerical scheme we have de- 
signed to solve the previous set of equations. We use the clas- 
sical "operator splitting" approach, decomposing the equa- 
tions into several steps that are solved in sequence. In the 
first step, called the "stellar source step" , we add all ionizing 
photons coming from stellar sources in the radiation field. 
In the second step, called the "transport step", we solve 
the hyperbolic system of conservation laws described in the 
last section. In the last step, called the "thermochemical 
step", we solve the right-hand side of our radiation trans- 
port model, together with the evolution of neutral hydrogen 
density and gas temperature. 



For the first step, we perform in each cell of the com- 
putational grid, indexed i, the following update: 

{N-,y^+^ ^{N^)^ + N;At (24) 

Here and in the foUowings, index n stands for the radiation 
field before the current step, and index n+1 for the radiation 
field after the current step. Since we have three intermediate 
steps, if we start a given time t with index n, we reach the 
next time step at time t + At with index n + 3. 



3.1 Transport Step 

The next operator we solve in our sequence is the hyperbolic 
system we have discussed in the next section: 



dt 



+ VF^ = 0, 



dt +''^^-'- 



(25) 
(26) 



In ATON, these equations can be solved either implicitly or 
explicitly, with the following integral form of the conserva- 
tion laws (expressed here in ID for sake of simplicity): 



At 



+ 



/2 



Aa; 

(-^7)1+1/2 ^ (-^7)1-1/2 



= 0, 



0. 



(27) 
(28) 



At Ax 

The flux function f = {Fj, P-,)'^ is evaluated at the intercell 
faces, indexed i + 1/2, and at time m, with m — n for an 
explicit scheme, whose stability condition imposes a strong 
constraint on the time step, or with m = n + 1 for an im- 
plicit scheme, unconditionnaly stable. Following 'Gonzale^ 
|et al.| ( |2007[ ), we compute the intercell flux using standard 
methods designed for computational fluid dynamics such as 
the Godunov method or, more generally, such as the class 
of upwind schemes. If we note W = (A''.^, ^^7)"^ the vector of 
state variable, the intercell flux depends on the left and right 
states with respect to the interface: 



1.- 



1/2 



■ ^ (Wi , I4i^ 



(29) 



Gonzalez et al. ( 2007 1 have tested various flux functions with 



respect to the Ml model, and came up with two possibilities, 
namely the Harten-Lax-van Leer (HLL) flux function, for 
which we have: 



)i + l/2 



\+ - A- 



(30) 



where A"*" = max(0,Amaa;) and A^ = min(0,Ami„) are the 
maximum and minimum eigenvalues of the Jacobian matrix 
of the system evaluated at the i-th and i -\- 1-th cells (for 
more details, the reader is encouraged to read |Gonzalez et al.| 
(2007[)) and the Global Lax Friedrich (GLF) flux function 



for which the maximum wave speed is taken equal to the 
speed of light: 



(31) 



As we will demonstrate in the test section, the GLF flux, 
more diffusive by nature, turns out to give results very simi- 
lar to the HLL flux, while being much simpler to implement, 
since it does not require to determine the eigenvalues of a 
rather complex hyperbolic system. 
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3.2 Implicit versus explicit ? 



out analytically: 



The time step of our scheme, At, is controlled by the 
Courant condition, which writes 3cAt/Ax < 1, if the in- 
tegration is performed explicitly (m — n). In cosmological 
problems, this can result into very small time steps, although 
the overall solution evolves quite slowly. The standard solu- 
tion is to use larger time steps, related for example to the 
much lower ionization front typical propagation speed, but 
we must rely on implicit integration, for which m — n + 1. 
In the latter case, we fix the Eddington factor x s^nd the 
various eigenvalues A to their initial value (index n), so that 
the implicit form of our solver becomes a linear system, as 
can be seen from Equations ( |21[ |, ( |27[ ) and ( |28[ ). This linear 
system is then solved using standard sparse solvers, among 
which the Gauss-Seidel solver gave the best results. 

Another alternative to overcome the limitation of the 
explicit scheme, while avoiding the complexity of the im- 



plicit scheme, was proposed by Gnedin & Abel (20011: the 



"Reduced Speed-of-Light Approximation". The idea is to 
replace the actual speed-of-light c by an effective speed- 
of-light c <C c in all the previous equations (including the 
photoionization rates), where 10~^ < c/c < 1 is typically 
considered. By reducing the speed-of-light, the Courant con- 
dition becomes looser and larger time steps can be consid- 
ered. Even though it introduces an approximation, no con- 
sequences is found as long as phenomena where the actual 
speed of propagation remains smaller than c are considered, 
like the expansion of an HII region. As shown in the follow- 
ing, using c instead of c has no consequence on the accuracy 
of the results. Furthermore, there are also a significant num- 
ber of cases where satisfying the strict Courant condition 
can be handled numerically, still avoiding the rather costly 
implicit solver. In this paper, we have never used the im- 
pUcit scheme, although in other physical conditions, using 
the implicit solver can be unavoidable. 



3.3 Thermochemical step 

In our operator splitting approach, the last step solves for 
the chemical evolution of Hydrogen and for the coupling be- 
tween radiation and matter. Physical quantities such as the 
gas temperature or the ionization fraction must be updated. 
The equations we solve here are the foUowings: 



dt 

i)r 

driHo 
dt 
de 
dt 



= —nuoCa^Nj + nonH^(aA — Q-b) (32) 

= — n-HoCa^F^ (33) 

= —nuoCa^Nj + ncUji^aA — n^nYiaP, (34) 

= n-agC(T^€^N^ — nHoWeAeHo — nn^rleAeH I . (35) 



In these equations, the coefficients a a, cxb, AeHq and AeH^_ 
all depend non-linearly on the gas temperature. The vari- 
ous chemical species' densities can be expressed as a func- 
tion of the ionization fraction x as = nu^. = xnu and 
riHo = (1 — x)nH- The gas internal energy is expressed as 
e — 3/2(1 + x)nnkBT . Because of the very small time scales 
involved, we solve this system using a fully implicit scheme. 
The first two equation are linear with respect to their main 
unknown, so that the implicit discretization can be worked 



n + l 



n+l 



1 + Atn'^+'- ea- 



rn 



(37) 



Injecting these equations into the rest of the system leads 
to an implicit system of 2 coupled non-linear equations to 
solve for the 2 variables x"'^^ and We describe in the 

Appendix the technical solution we propose to solve this 
problem. After this final step, we find the new thermochem- 
ical state (a;,T)"+^, and the corresponding new radiation 
state {N.y,Fry)"^^ , at the new time step t + At. 



4 TESTS AND RESULTS 

In their cosmological radiative transfer codes comparison 
( 2006 1 give a set of standard tests that 



Iliev et al 



project, 

were used to probe the validity of ATON and its implemen- 
tation. First the expansion of a Stromgren sphere is simu- 
lated, then the shadowing of an I-front by a dense clump is 
discussed, and a first attempt to model the propagation of 
radiations in a static cosmological field is presented in the 
end. 



4.1 HII region expansion with a constant 
temperature 

First, the classical situation where a single source emits a 
ionizing radiation that propagates through the surround- 
ing neutral medium is presented. It results in the classical 
picture where an ionized sphere , centered on the source, 
expands at given rate. As the recombination process starts 
to counter-balance the ionization, the front's position slows 
down and even stops, achieving a stationary regime. Assum- 
ing the source emits A''.^ ionizing photons per unit time the 
Stromgren radius, i.e. the stationary radius, is given by: 



rs 



3n; 



'inaBiT)njj 



1/3 



(38) 



where as(r) stands for the case B recombination rate at 
a temperature T and hh is the surrounding gas number 
density. The time evolution of the I-front 's position is given 
by 



riit) =rs(l- 



(39) 



where t,. is the characteristic recombination time given by 
1/tr — ctB{T)nH. From Eq. 39 one can see that the expan- 
sion slows down for t ^ tr. 

The setting for the numerical experiment is similar to 

The source is located 



al. 



(20061 



the one given by Iliev et 
at the corner of the simulated box and emits A'^^ = 5 x 
10** photons per second. The surrounding hydrogen has a 
number density of uh = 10~'^cm~^ with an initial ionized 
fraction x = 1.2 x 10~^ and its temperature remains fixed, 
with T — W^K . The simulation is performed on 64'' grid, 
with a 6.6 kpc box size. Refiexive boundary conditions are 
assumed. The source is switched on at t=0. The results are 
presented in Figs, [l] and [2] 
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Figure 1. Stromgren sphere test. GLF intercell flux has been used for the left panel's computations, HLL intercell flux has been used 
for the right panel's calculation. The lines stand for the ionized fraction profile and the neutral fraction profile, computed at t = 35 
Myr (top row) and t= 500 Myr (bottom row) with a reduced speed of light c = c/1000 (dot-dashed line), c = c/100 (dotted line) and 
c = c/100 (plain line). Curves are superimposed at t= 500 Myrs. 




Figure 2. Stromgren sphere test. The dotted (resp. plain) lines stand for the I-front position at different times (in units of the recom- 
bination time), computed with a reduced speed of light c = c/1000, c/100, c/10 (resp. dot-dashed, dotted and plain lines). HLL flux 
calculation is shown on the right panels, GLF flux calculation on the left panels. Bottom rows: ratio of the simulated front position to 
the analytic value. Top rows: I-front position evolution in units of the box size. The Stromgren radius is 5.4 kpc and the recombination 
time is 122.4 Myr. 



Fig.[T]shows the profiles of the ionized and neutral frac- 
tion computed with the two different types of intercell flux 
at t=30 Myr and t=500Myr. Clearly the two methods agree 
for an effective speed of ligiit that varies from c = 0.001c to 
c = c. For instance, both methods return an underestimated 
I-Front position as c is smafler than the physical propaga- 
tion speed of the front. Conversely, when c is large enough 
(typically c 0.01c), the two fluxes approximations agree 
and achieve convergence regarding the result. The same cal- 
culation is performed in |Iliev et al.| ( |2006[ ) and the results 
shown in Fig. [T] are consistent with the profiles found by the 
codes that took part to their comparison project. 

In Fig. (|2|, the I-front 's position is compared to the the- 
oretical calculation. From the top panels of Fig. H, it clearly 



appears that the current code accurately reproduces the fast 
expansion of the ionized bubble at early times and the slow 
down as the I-front position gets closer to the Stromgren ra- 
dius. At later times, a stationary state is achieved and the I- 
front stops. Again, the GLF and HLL calculation agrees and 
return a final radius 4% larger than the theoretical Strom- 
gren radius. The comparison to |IIiev et al.| ( [2006) shows that 
it is consistent with the other types of calculations which 
tends to overestimate the final radius' value. The bottom 
panels show the comparison between the simulated I-Front 
position and the theoretical one at each time step. Clearly 
the two types of intercell fluxes (GLF and HLL) return sim- 
ilar results and both present a drift toward greater values of 
the front position. For t=500 Myr, the discrepancy is 4%. 
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Interestingly, this drift is similar for all values of c and the 
results only differ by the time required for the profile to 
'converge' toward the drift. For instance, the c = 0.001c cal- 
culation requires ~ tree to catch up the main drift, while it 
is instantaneous for c = 0.1c. The results found in llliev et al.l 



1 2006 1 shows a similar behavior. 



4.2 HII region expansion with a variable 
temperature 

In the second test, the HII region expansion is investigated 
while the gas temperature is allowed to vary as a conse- 
quence of photoheating. The numerical setup is identical to 
the previous test with two exceptions. First, the source is 
now a 10^ K black body source while the flux of ionizing 
photons remains unchanged. Second, the initial gas temper- 
ature is set to 100 K and its evolution is modeled according 
to the Eq. ( |19[ ). The source is turned on at f = and the 
simulation is run over one recombination time with c = 0.1c. 

Fig. ([3| presents the temperature and neutral fraction 
maps after 100 Myrs and Fig. Q presents the time evolu- 
tion of the temperature and ionized fraction profile. A quick 
comparison between the constant temperature experiment 
and the current one indicates that the results qualitatively 
agrees. For instance, taking the t=35 Myrs profile as a ref- 
erence, the transition region is in both case located at a 
radius equals to 0.45 , in units of the box size. However the 
transition operates on a larger range of radii compared to 
the isothermal case, the background neutral fraction being 
achieved at r=0.8 box length at t=35 Myrs while this re- 
gion does not extend further than r=0.65 in the constant 
temperature case. This relates to the non-monochromatic 
source used in the current case : the average photon-energy 
is larger than previously, hence the average cross-section is 
diminished and this leads to a smoother transition from the 
fully ionized to the neutral regime. The temperature pro- 
file shows a steady evolution of the limit between the hot 
medium (TJs; lOOOOK) and the cold neutral medium. The 
transition between these two regimes is much sharper than 
observed for the neutral fraction profile. Furthermore the 
inner temperature profile is practically fiat (except close to 
the center) while the neutral fraction values spans over a few 
orders of magnitude. It can also be observed on the maps 
where no dynamics can be seen on the temperature maps 
while a gradient in the neutral fraction is clearly visible. 

Comparisons with Iliev et al. ( 2006 1 results show that 



ATON's calculations are consistent with other codes, how- 
ever the transition region's profile between ionized and neu- 
tral gas is sharper in the current calculations than returned 
by the other codes. It could be related to the current spec- 
tral treatment, where a single photon-population (having 
an energy representative of the source spectrum) is taken in 
account, the transition's profile refiecting its typical energy. 
Conversely, a more complete treatment of spectral harden- 
ing is expected to produce smoother profiles, as high-energy 
photons have larger mean free-path and travel further in the 
neutral hydrogen. 

4.3 Shadowing by a dense clump 

The third test investigates the code's ability to deal with 
density clumps along the photon's path. Such clumps are 



expected to slow down the I-fronts propagation and cre- 
ate regions which are shielded from the radiations due to 
the shadow trailing behind high density regions. Such situ- 
ations are likely to happen in a cosmological context where 
shielding is provided by e.g. gaseous structures or filaments. 



The setting suggested by Iliev et al. (20061 is adopted: a 



6.6 kpc box is considered with a dense and cold clump on 
the radiation's path. The background gas density is p = 200 
m"'' with a initial temperature T = 8000 K. The clump is 
sphere located at (5, 3.3, 3.3) kpc with a radius r — 0.8 kpc. 
Its density is 40000 m~^ and its temperature is T = 80 K. 
A stationary photon fiux ^ — 10^''m~^s~^ photons is ig- 
nited at t = with a typical energy corresponding to a 10^ 
K black-body source. Calculations were performed with the 
HLL intercell fiux with c = c. 

Neutral fraction and temperature maps and profiles are 
given in Fig. ([7| and Fig. (|5|. Clearly shadows are created in 
the clump trail: the gas remains neutral behind the clump 
and its temperature remains lower than the surroundings 
due to the lack of heating photons in this region. The neu- 
tral fraction is almost zero outside the clump due to the 
incoming flux of photons ad rises to a constant level around 
X ~ 0.01 on the 'enlightened' side of the clump as light 
penetrates the over-density. As illustrated by Fig. Q, the 
I-front propagates deeper into the clump and after a fast 
propagation in the optically thin medium, the front almost 
stops and progress slowly trough the dense medium. On the 
other hand, the temperature proflles shows a similar behav- 
ior with a sharp transition region from ~ 10500 — 11000 K 
to 80 K, the initial clump's temperature. In the 'trailing' 
region, temperature is initially lower than in the leading re- 
gion. However, it is clear from the maps in Fig. ([7| that the 
limits of the trailing shadow are not parallel to the direction 
of the incident flux, and a difi'usion cone appears. As time 
advances, the temperature in the initially shielded region 
rises, until it reaches a level close to the one observed in the 
exposed gas. Finally, Fig. ([8| describes the time evolution 
of the average ionized fraction and temperature inside the 
clump. Both curves shows a fast increase during the first mil- 
lion years followed by a shallower evolution toward x = 0.8 
and T = 11000 K. Up to t=10 Myrs, these evolutions are in 
a quantitative agreement to the ones presented in |Iliev et al.| 
( 2006]) , which indicates that the current scheme captures the 
overall picture of the I-front trapping, even though the lack 
of multi-frequency treatment leads to some discrepancies in 
the detailed description of the process. 

From a physical point-of-view, a certain level of diffusiv- 
ity is expected : atoms recombine and radiate in an isotropic 
fashion. In the current scheme, isotropy is taken in account 
by the spherical component of the Eddington's tensor and 
its relative contribution is set by the ratio of the local flux by 
the energy ( the reduced fiux /). From the maps in Fig. 
the clear conic shape of the trailing shadow appears as a 
manifestation of this isotropy. However, ATON's predicts 
a larger shadow's extent compared to experiments in |Iliev| 
et al. ( 2006 1 where it remains parallel and clear cut : this 



discrepancy is likely to be related to the lack of modelisa- 
tion of the recombination's isotropy. As stated in the first 
section, ATON's can mimic the behavior of a ray-tracing 
code by adding an extra artificial term in the flux equa- 
tion (see Eq. Q and Eq. (??))• The computations using 
this scheme are significantly different from the one using 
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Figure 3. Stromgren sphere test with non uniform temperature. Left : Neutral fraction map in the source plane. Right ; temperature (in 
Kelvin) map in the source plane. These maps were computed 100 Myr after the central source (located here in the bottom left corner) 
has been switched on. The experimental setting is similar to the constant temperature test, described in sect ion [4. 1[ 




Figure 4. Stromgren sphere test with non uniform temperature. The left panel shows the average temperature profile while the right 
panel shows the the average ionization and neutral fraction profile. These profiles were computed lO(plain), 35 (dotted) and 100 (dashed- 
dotted) Myrs after the central source has been switched on. The experimental setting is similar to the constant temperature test, described 
in section 



the regular scheme and are in good agreement with results 
shown in Iliev et al. (20061 . In Fig. ([7|, it clearly appears 



that the 'ray-tracing' scheme is much more 'conservative' in 
terms of flux geometry: the shadows are clear cut behind the 
clump. It is particularly striking in the temperature maps 
where basically no shadow persists after a given duration in 
the regular calculation while it remains and exhibits sharp 
edges in the current enhanced calculation. The same effect 
can be noted in the neutral fraction maps where the cone- 
shaped shadow is replaced by a straight cylindrical one. It 
persists over the 15 Myrs of the calculation. In terms of pro- 
files, it can be seen from Fig. Q, Fig. ([8| and Fig. (|6| that 
the I-front propagation is affected too and their progression 
through the clump is faster. The radiation's flux is more 
'rigid' and the influence of diffusion is lowered : the ion- 
ization of a propagating front is more efficient especially in 



dense regions where I-front can be significantly slowed down. 
Also, the smaller diffusivity keeps the initial temperature 
profile in the shielded region unchanged. Let us emphasize 
that these results are obtained using an artificial flux term 
and should be considered with caution since diffusivity must 
appear at some point. As shown hereafter, the differences re- 
main small in realistic cosmological situations but still : in 
the current setting clear cut shadows are obtained by sup- 
pressing the isotropic recombination and as shown in the 
results obtained using the regular ATON's scheme, results 
can be significantly different in certain situations. 

The other differences with the calculations presented in 



Iliev et al. ( 2006 1 can be explained by the lack of multi- 



frequency treatment. Because the scheme do not take in ac- 
count the effect of high frequency photons with large mean 
free-path, no preheating on large distance is being modeled. 
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With the current single-energy treatment, the effect of these 
photons tend to be underestimated and all processes (ion- 
ization and photo-heating) occur on small scales, leading to 
the observed sharp profiles. In Fig. ([5| , the ionization frac- 
tion transition from ionized background to neutral clump is 
much sharper in ATON's case and no steady decrease can 
be observed. Other codes do not present a sharp drop in ion- 
ized fraction and a smoother transition exists between the 
neutral fraction of the irradiated side to the shielded side. A 
similar discrepancy exists in the temperature profile, where 
the drop in temperature occurs on smaller distances. The 
average ionized fraction within the clump is also a few per- 
cent smaller than found in other calculations at late times. 
In other words, the progression of the I-front within the 
clump appears to be slightly slower in ATON's calculations, 
the trapping is more efficient. It results from the lack of 
preheating and ionization by high energy photons. In the 
prospect of cosmological studies, where sources are likely to 
be embedded in high density regions, the radiation is there- 
fore expected to take a slightly larger amount of time to 
escape from the haloes in the current scheme compared to 
other calculations. 



4.4 Static Cosmological Field 

The final test consists in modeling the propagation of I- 
fronts in an heterogeneous medium with multiple sources. 
The gas density field consists in a single snapshot extracted 
from a 128^ cosmological simulation. The comoving box 
length is 500ft~^ kpc and for sake of simplicity only the 
z = 8.9 snapshot is taken in account and the medium is 
considered as being static, both from the point of view of 
cosmology (e.g. no overall expansion) or local motion. A sim- 
ple source model and distribution was also derived, based 
on the properties of the biggest haloes. All the sources are 
turned on a the same time and have lO^K black-body spec- 
trum. The initial temperature is 100 K and the simulation is 
ran over 400 000 years, i.e. well before the stationary regime 
achieved. ATON's computations are made with an effective 
speed of light equal to c. Further details can be found in 



Iliev et al |(2006| ). 

Fig. (|9 1 and Fig. ( 10 1 present the neutral fraction and 
the temperature map at t=0.05 Myr and t=0.4 Myr : only 
the mid-section planes of the simulation are shown. These 
maps were obtained with the two versions of ATON (GLF 



and HLL) and compared to the same tests performed by 


FTTE (Razoumov & Cardall 


(20051) and C2RAY CMellema 


et al. (2006l) as part of the 


Iliev et al.|(|2006) comparison 



project. Let us emphasize that the simulation is run over 
400 000 years only, i.e 0.1% of the recombination time and 
the regime investigated in this test is highly non-stationary. 
Clear differences can be noted. First, ATON's calculations 
seem to present a delay in the ionization propagation. It can 
be seen from both maps where ionized and hot regions are 
typically less extended around the sources in ATON's cal- 
culations compared to the other results. Second, hot regions 
are more tightly correlated to ionized regions in ATON's 
scheme, while FTTE and C2ray show extended and com- 
plex structures of hot gas. Since the transition between hot 
and cooler regions is more abrupt in the ATON's calcula- 
tions, it relates to the more compact ionized regions that the 
current scheme predicts. Furthermore, let us recall that the 
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Figure 5. Shadowing by a dense clump using HLL intercell flux. 
Top: Temperature profile along the middle section of the box. 
Bottom: Neutral and ionized fraction profile along the same sec- 
tion. Profiles were computed at t= 1 (bottom row), 3 (middle 
row) and 10 (top row) Myrs. The radiation comes from the left 
side of the maps. Thin solid lines stand for the computation us- 
ing the 'regular' Ml scheme while thick dashed lines stand for the 
computation using our "ray-tracing" scheme. 



current implementation of the scheme is "monochromatic", 
even though the average photon's energy refiects the spec- 
trum of the sources. As a consequence, fronts are sharper 
and high-energy photons with large travel distances are not 
available and cannot preheat and preionize the gas on large 
scales. 

Still, given this strong restriction to a single frequency- 
group treatment, results can be seen as qualitatively satisfy- 



ing. The I-front positions (see Fig. (Ill) with ATON share 
similar features with the other codes. For instance they co- 
incides along high density regions such as filaments axes. 
Densities are such that the frequency dependence of cross- 
section do not penalize the scheme at the current resolution 
compared to the others and conversely, I-fronts in voids are 
"late" in ATON: this supports the fact that the differences 
are due to the current simplified spectral treatment. This 
is particularly evident in the 0.05 Myr map, where a good 
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Figure 6. Shadowing by a dense clump. The I-front position as a 
function of time for the regular ATON's scheme (thin solid line) 
and "ray-tracing" scheme (thick dashed line). The position of 
the I-front, relative to the clump center is given in units of the 
Stromgren length defined as £g = F/p^ag. 



agreement coincides with filaments and a poor agreement 
coincides with voids. The 0.4 Myr maps still hint such cor- 
relations but they appear as less obvious since it results from 
a longer evolution. 

The overall ionization evolution is presented in 
Fig. (12 1. The average ionized fraction is defined by 



N ' 



(40) 



where A*' stands for the total number of grid cells, i is the 
cell index and the mass- weighted average ionized fraction by 



(41) 



where pi is the local gas density within cell i. The ratio of 
these two quantities is equal to the ratio between the average 



gas density of ionized regions to the average density (Iliev 



et al. (20061) 



(42) 



-^^ionized 

onizcd 

A ratio larger than one implies that ionized regions have 
densities larger than the average. Conversely a ratio smaller 
than unity means that voids dominate the 'population' of 
ionized regions. Compared to the other codes calculations, 
ATON exhibits an overall shift of Xm{t) and Xv{t) in am- 
plitude (~ 4% with FTTE and ~ 12% with C2RAY) while 
the global trends are similar. One can notice that for all the 
codes the curves for Xm and x^ cross each other at some 
point, implying that high density regions are ionized first 
while voids tend to be ionized later. The switch operates 
later in ATON and confirms the code's slowness in low den- 
sity regions. Still the discrepancy with other codes remains 
limited and is encouraging in the prospect of a more com- 
plete treatment of multi- frequency transfer. Let us also em- 
phasize that the current set up investigate a highly transient 
regime and the previous tests have shown that a much better 
agreement is expected on longer time scales. 




CS o 




Figure 7. Shadowing by a dense clump using HLL intercell flux. 
Top: neutral fraction (top row) and temperature map (bottom 
row) along the middle section of the box, measured at t=lMyr. 
Bottom: same measurements at t= 3 Myrs. Left columns stand for 
computations with a "ray-tracing" scheme, while right columns 
stand for the 'regular' Ml scheme. 



Finally, all these tests were performed using also our 
"ray-tracing" scheme and no significant differences were 



found (see Fig. ( 13 \). The averaged ionized fractions remain 
unchanged at the percent level. The maps and the i-front po- 
sitions (not shown here) are quasi-identical. It implies that 
even though the isotropic recombination is artificially sup- 
pressed, it does not have a strong impact compared to the 
lack of multi frequency treatment for large scale calculations. 
In the context of cosmological reionization, future develop- 
ments would nevertheless focus on these aspects first. 



5 DISCUSSION & CONCLUSION 

ATON is a Eulerian scheme for radiative transfer that re- 
lies on a momentum description of the radiative transfer 
equation. The hierarchy set up by the conservation equa- 
tions of energy and radiative flux is closed by means of a 
relation between the radiative pressure and the energy. The 
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t [Myr] 

Figure 8. Shadowing by a dense clump. The time evolution of 
the average ionized fraction (top) and the average temperature 
(bottom) within the clump for the regular ATON's scheme (thin 
solid line) and the "ray— tracing" scheme (thick dashed line). 



'Ml' closure relation has been retained : it expresses the 
Eddington tensor as the combination of a pure transport 
configuration of the radiation and a purely diffusive geome- 
try. In the intermediate cases , the relative contribution of 
these two regimes is obtained from the local flux properties. 
This scheme is complemented with a a simple treatment of 
the local chemistry, to compute the neutral hydrogen abun- 
dance, and the conservation of energy, to compute the pho- 
toheating/cooling. At the current stage, the multi-frequency 
treatment has not been implemented, even though there are 
no a priori difficulties for such an extension. ATON is a 
simple scheme that will eventually be used to investigate 
astrophysical ionization processes such as the cosmological 
reionization. 

ATON has been tested following the experiments sug- 



gested by Iliev et al. (2006j): the propagation of an HII re- 
gion, with and without auto-consistent heating, the shadow- 
ing by a dense clump and the I-front propagation in a static 
cosmological field. Most of the results obtained by ATON are 
in agreement with the calculations presented in this article. 
The main differences arises from the lack of multi-frequency 
treatment : because the spectrum hardening is not taken in 
account, ionization and heating processes occurs on small 
scales and for instance no large distance pre-heating due to 
high energy photons is being modeled at the moment. It 
results on a loss of radiative energy on small scales and I- 
fronts tend to be slower in low density regions. However, it 
only causes problems in highly non stationary regimes while 
ATON catches all the details of the processes in situations 
where I-front propagation are slowed down. The issue of the 
suppression of isotropic recombination was also assessed. For 
this purpose, we used a modified scheme that mimics ray- 
tracing codes by taking in account an anisotropic source of 
flux due to recombination : the results obtained are very 
similar to ray-tracing or Monte-Carlo codes. However, little 
differences were observed in the cosmological test, empha- 
sizing the greater impact of the mono-group treatment. 

Among the current and future developments, the multi- 





Figure 9. Static cosmological field test. Maps of the neutral frac- 
tion computed at 0.05 Myr {Top) and 0.4 Myr {bottom) after 
the sources were ignited, the four calculations were made using 
ATON-GLF, ATON-HLL, C2RAY and FTTE. 



frequency treatment will follow in order to investigate accu- 
rately ionization processes in highly non stationary regimes. 
Such situations would occur on small scales close to the 
sources, and these studies would be valuable to constrain 
e.g. the escape fraction on scales smaller than the resolution 
of large volume simulations. Also, ATON is limited to the 
post-processing of simulations and lacks the coupling that 
exist between the dynamic of the gas and the radiation. Be- 
cause of its Eulerian nature, the current scheme can be eas- 
ily coupled with grid based hydrodynamical codes and it is 
planned to bo part of the AMR cosmological code RAMSES 
(Jteyssier 2002). 

Among the astrophysical applications of ATON, the 
study of the cosmological reinsertion is a primary objective. 
The post-processing of the AMR and SPH version of a 
large ongoing hydrodynamical simulation is on the way and 
will allow to compare the impact of the source distribution 
and the gas geometry on the reionization process generated 
by each version . It would also lead to predictions on 
the geometry of the 21 cm emission, in the prospects of 
experiments such as SKA or LOFAR. Comparisons on 
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Figure 10. Static cosmological field test. Maps of the temper- 
ature computed at 0.05 Myr {Top) and 0.4 Myr (bottom) after 
the sources were ignited, the four calculations were made using 
ATON-GLF, ATON-HLL, C2RAY and FTTE. 



Figure 11. Static cosmological field test. I-front positions (de- 
fined as a 50% ionized fraction) computed at 0.05 Myr (Top) and 
0.4 Myr {bottom) after the sources were ignited, the four calcu- 
lations were made using ATON-GLF, ATON-HLL, C2RAY and 
FTTE. 



the propagation of radiation in and around the biggest 
objects, at high resolution are also on the way in order to 
investigate the impact of small scales clustering of the gas 
and the star distribution. Hopefully, with the improvements 
mentioned previously, the current scheme will be able to 
investigate all the astrophysical processes where radiation is 
relevant, while remaining simple to conceive and implement. 
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Figure 12. Static cosmological field test. Top : the evolution of 
the average ionized fraction (thin lines) and mass-weighted av- 
erage Xm (thick lines) . These quantities were computed on snap- 
shots produced by C2RAY, FTTE, ATON-GLF and ATON-HLL. 
Bottom : the evolution of the ratio Xm/xy for the same four cal- 
culations. 



Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, 

New Astronomy, 11, 374 
Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006, 

MNRAS , 372, 679 
Razoumov A. O., CardaU C. Y., 2005, MNRAS , 362, 1413 
Teyssier R., 2002, AAP , 385, 337 

Toro E., 1999, Riemann Solvers and Numerical Methods for 
Fluid Dynamics. Berlin, Germany, Springer- Verlag, 1999, 
624 p. 



APPENDIX A: IMPLICIT COMPUTATION OF 
IONIZATION FRACTION 

In ATON, the coupling between radiative energy and ioniza- 
tion fraction is solved implicitly. Let us call A'^^ the density 
number of photons (i.e. the energy density in single photon 
units), F-y the associated flux, P-y the radiative pressure and 
X the ionization fraction. Let us also define p as the ini- 
tial hydrogen density and a, as and /3 as respectively the 
case A and case B recombination rates and the collisional 
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Figure 13. Static cosmological field test. The evolution of the 
average ionized fraction x„ (thin lines) and mass- weighted average 
Xm (thick lines) using the regular scheme and the "ray— tracing" 
scheme (with flux mention). An offset of 0.02 has been applied 
to the latter quantity. GLF (top) and HLL (bottom) intercell 
flux were used. Clearly the "ray-tracing" scheme do not lead to 
significant changes in the global evolution of the ionization. 

ionization rate. First we consider a fixed temperature (from 
the previous timestep) for the rates computation. The set of 
coupled (ID) equations is given by 

AN-, AF^ 
~AF ~d7 



AF-, 



AP^ 



At ~^ Ar 



N; + N'^'"' ~ pa-.,cN.y{l- x) (Al) 
-pa-ycFj. (A2) 



Here A''* and A'^^'^'^ stands for the point-like (namely the 
'stars') and the diffuse source of ionizmg photons due to 
recombination and stands for the photo-ionization cross- 
section averaged over the A^.^ spectrum. The ionization equa- 
tion is given by 

= ap'a;'^ - Px{l - x)p^ - p{l - x)N^a^c. (A3) 

The equation over n can be rewritten as: 

-^ + -^ = S-aBpx^+Px{l-x)p^-p—, (A4) 

where we assumed that oip^x'^ — N^'^'^ + otBp^x'^ . With p 
labeling the time step index, we write ANj/At = (A^^^^ — 
NFf)/At and x = x'' and X — x''^^ and the implicit formu- 
lation of Eq. ||A4l) is given by: 



A'^+^ = N^' + /3p^(l - X)XAt - agp^X^At - p[X - x), (A5) 
where TV' is the explicit solution of the pure advection equa- 



tion (i.e. Eq. (All with a zero r.h.s.) given by: 



A^; ^A-^At- ^At-f A-^. 



(A6) 



Combining Eq. (All and Eq. |A5l, one can obtain a third- 



order polynomial Q^{X) = ■mX^~+ nX^ +pX + q, where the 
coefBcients are given by 

m = (aB+ l3)p^At 
n — p 



a + l3)p/a^c - asp At - 2(3 p At 



(A7) 
(A8) 
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p = -p{l + x) - N!y - l/a.ycAt + I3p/a.yc + l3p^At{A9) 
q = N!y + px + x/a-,cAt, (AlO) 

and where X is given by finding the roots of Q3{X). Know- 



ing X, N^, is obtained from Eq. ( A5I and F-y can be derived 
from 

F' 

^^^^ " 1 + ca^pAt{l - X) ' ^^^^^ 

where Fl^ is the explicit solution of the the pure advection 
equation of the flux, given by: 

dP^ 
dr 

The temperature at time t + At is obtained a posteri- 
ori, having computed values of A'^^"'"^, FJ^'^^ and x^^^ while 
using the temperature set at the previous time step. From 



Eq. ( 19 1, the (discrete) equation that rules the temperature's 
evolution is given by; 

rpp+i _j'P 2{H-C - |p(l + X)kBTP+\X - x)/At) 



At 3p(l + X)fcfl 

Ti and C being the heating and cooling rates. Therefore, 



(A13) 



Eq. (A13l can be solved separately in an explicit manner 
where the R.H.S depends on instead of T^+^ . This equa- 
tion is purely local and is solved at each cell's location with- 
out any spatial coupling. By simple algebraic manipulation, 
the updated value of the temperature T^+^ can in principle 
be easily obtained. However, the cooling time becomes much 
shorter than At as temperature reach typical value of ~ 10'* 



K, therefore Eq. (A13l cannot be solved without resolving 



this typical time scale. 



In practice, Eq. (A13l is sub-cycled during a radiative 
time step, using At = 0.9tcooi- At each 'temperature' itera- 
tion, the cooling rate is updated, setting a new At for the 
next iteration. As a consequence such a procedure can sub- 
stantially increase the computing time. This increase can 
be reduced by stopping the sub-cycling when some relative 
convergence of the temperature is achieved. In the tests pre- 
sented hereafter, a condition such as a 10~® convergence af- 
ter at least 100 sub-cycles is found to give satisfying result in 
terms of ionization fraction distributions. Finally the whole 
process is repeated: the ionized fraction is re-computed us- 
ing the updated temperature until a global convergence is 
achieved. 

This model is admittedly oversimplified and a full im- 
plicit treatment of ionization and photoheating would be 
preferable. However, it appears from subsequent experi- 
ments that taking in account the complexity of the energy- 
ionization coupling with a temperature that varies in back- 
ground returns satisfying results. Since temperature is es- 
sentially crucial to compute recombination rates, which do 
not depend strongly on temperature in our case, this simple 
model appears to be accurate enough at the current stage. 



